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We revise the Thomas-Fermi approximation for describing vortex states in Bose condensates of 
magnetically trapped atoms. Our approach is based on considering the ft ^ limit rather than the 
N ^ oo limit as Thomas- Fermi approximation in close analogy with the Fermi systems. Even for 
relatively small numbers of trapped particles we find good agreement between Gross-Pitaevskii and 
Thomas-Fermi calculations for the different contributions to the total energy of the atoms in the 
condensate. We also discuss the application of our approach to the description of vortex states in 
superfluid fermionic systems in the Ginzburg-Landau regime. 
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I. INTRODUCTION 

The discovery of Bose-Einstein condensation in 
trapped alkali-metal gases at ultra-low temperature [1-^3] 
has developed a huge amount of experimental and the- 
oretical investigations. The experimental conditions are 
such that the atomic gas is at very low density and that 
the interactions can be parametrized in terms of a scat- 
tering length a. In this situation a mean-field descrip- 
tion through the Gross-Pitaevskii equation (GPE) [4, 5] 
is able to give, at least at low temperature, a precise de- 
scription of the atomic condensates and their dynamics 
[6-9]. 

One important question concerns the superfluid char- 
acter of the Bose condensates. Among other properties, 
the existence of quantum vortices is a signal of the su- 
perfluidity. The possibility of trapped quantized vortices 
was one of the primary motivations of the GP theory 
[4, 5] and some amount of theoretical work about this 
topic has been developed during the last years [6, 7, 9- 
12]. The experimental evidence of such quantized vor- 
tices has recently been verified [13, 14]. 

Since the number N of atoms involved in the con- 
densate is generally large, it is natural to think that 
the Thomas-Fermi (TF) approach can be applied exten- 
sively in some aspects of the Bose-Einstein condensation 
in traps. This TF limit is usually identified with the limit 
of number of atoms N going to infinity rather than to be 
interpreted as the ?i — s- limit as it happens in the case 
of Fermi statistics. Recently, the TF approximation for 
the ground state of Bose-Einstein condensates of mag- 
netically trapped atoms has been discussed as the ?i — > 
limit [15]. From this point of view the TF kinetic energy, 
which is dropped in the N ^ oo limit of the ground-state 
calculation, can be obtained for any number of particles. 
In this h ^ limit, a good agreement between the GP 
and TF kinetic energies is found even for low and inter- 
mediate number of particles. With the interpretation of 



the TF approach as the h ^ limit, it is also possible 
to perform semiclassical TF calculations for the ground 
state of Bose-Einstein condensates of atoms with nega- 
tive scattering length (^Li atoms) and to compute the 
excitation energy of collective monopole and quadrupole 
oscillations where the kinetic energy of the ground state 
of the condensate plays a crucial role [15]. 

The TF limit considered as ^ oo limit has also been 
applied to the description of vortex states [10, 11, 16, 17]. 
In this case it is assumed that the radial and axial kinetic 
energies can be neglected, and only the rotational kinetic 
energy is retained [16]. This approximation, however, 
gives a bad description of the vortex-core region. A better 
description of this region in the limit of large A^ can be 
achieved by splitting the condensate wave function into a 
product of a slowly-varying envelope, which is obtained 
by completely neglecting the kinetic energy, times the 
solution of the GPE describing a vortex in homogenous 
matter [18]. In contrast to these large- A^ methods, in 
this paper, we will again consider the TF approximation 
as the h ^ limit in order to describe the vortex state 
semiclassically. In addition to the formal aspects, this 
approach has the practical advantages mentioned above, 
i.e., that one can calculate the kinetic energy, that it can 
therefore be used also in the attractive case, and that it 
works well also in the case of relatively few particles. 

The experimental and theoretical achievements in 
Bose-Einstein condensation have also triggered the inves- 
tigation of trapped Fermi gases at very low temperatures 
[19-21]. One of the most important goals of the exper- 
iments is to reach the BCS transition to the superfluid 
phase, associated with the appearance of a macroscopic 
order parameter of strongly correlated Cooper pairs in 
dilute gases of trapped fermionic atoms. Several theo- 
retical studies about this topic have recently been devel- 
oped [22, 23]. In the case where the critical tempera- 
ture is much higher than the spacing between the lev- 
els in the trap, the macroscopic order parameter can be 
obtained through the Ginzburg-Landau equation (GLE) 
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[24], which is formally equivalent to the GPE. In a recent 
publication [25] also vortex states were discussed within 
the framework of the GLE. As a second application of 
our TF approach we will briefly discuss vortex states in 
a superfluid gas of trapped fermionic atoms. Due to the 
analogy between the GPE and the GLE our semiclassical 
approach can immediately be transferred to this problem. 

The paper is organized as follows: In the second sec- 
tion wc establish the TF theory projected on states of 
defined z component of the angular momentum and ap- 
ply it to describe vortex states of a non-interacting Bose 
condensate. In the third section wc include the inter- 
action among the atoms in the trap and compare our 
semiclassical prediction with the results obtained from 
the quantal solution of the GPE for several typical ex- 
amples. The fourth section is devoted to the discussion of 
vortex states in superfluid trapped Fermi systems. Our 
conclusions are laid out in the last section. 



II. THE THOMAS-FERMI APPROXIMATION 
TO STATIC VORTEX STATES 

We start by considering states having a vortex line 
along the z axis and all the atoms flowing around it with 
quantized circulation. The order parameter can be writ- 
ten in the form [6, 11] 



$(r) = <Afe,^)e-*', 



(1) 



where and z are the radial and axial coordinates, 

is the angle around the z-axis, k is an integer, and 
(^(r^ , z) = p(r^ , z), p{r^ , z) being the density. The vor- 
tex state has a tangential velocity v = hK/{m,r^) where 
K is the quantum of circulation, and the angular mo- 
mentum along the z axis is Nhn. The function ^(r^ , z) 
is obtained as the solution of the following non-linear 
Schrodinger equation 
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2m\dr^ ^1 dz"^/ ' 2mrf 
+ Vecot{rj_,z)+g(f)'^{r^,z) (l>{r^ , z) = n (l>{r^ , z) , (2) 



which is the GPE for the static vortex state problem. In 
Eq. (2), Vext is an external potential which for simplicity 
we have considered to be a spherical harmonic oscillator 
(HO) with frequency co, 



VeA^ = ]^muj\rl+z^). 



(3) 



The coupling constant is given hj g = Anti^a/m with m 
the atomic mass and a the s-wave scattering length. 

For the remaining part of this section, we will concen- 
trate on the non-interacting case, i.e., V{r) = Vext{r)- 
The effect of interactions will be considered in the next 
section. For non-interacting particles one recovers the 



case of a stationary Schrodinger equation for the har- 
monic oscillator potential, which is solved by 



(4) 



with the HO length uho defined by qho = \/h/{mu)). 
The corresponding energy eigenvalue is given by 



(5) 



To derive the TF approximation to the quantal solu- 
tion of the non-interacting vortex state (4), we start from 
the complete set of eigenfunctions oip^, Pz and 

(rlfc„fc^,«) = J,(A;,rJe^'"^e''=^^ (6) 

which are normalized to 
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d^r {r\k„k^,K) {k'^,k[,K'\r) 



= S{k^-k'J6{k,-k',)d,,,. (7) 

At lowest order in 7i (i.e., at TF level), the corresponding 
single-particle propagator [26] can be written as 



where R = (r -I- r')/2. Eq. (8) has been obtained under 
the assumption that all the gradients of the potential can 
be neglected, which is the usual hypothesis of the TF 
theory. From now on we restrict ourselves to some given 
value of K. The spectral density matrix is easily obtained 
as the inverse Laplace transform of the propagator [26]: 
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k^Mk^rJMky^e'^'^-'^'^ 



X e'"' (^--') 6(^^^-V{R)- ^%^) • (9) 

Its local part, g^{f) = g^{f,f), is proportional to the 
density of the Bose condensate. After performing the fc^ 
integral we obtain for the density 
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dk. J2 [^k^oir^-kl r^) 0[fi - V{r)] , 



(10) 
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where 



v/2m[M - V{r)] 



(11) 



and is the normalization constant. 

As it is done for the ground state [15], is determined 
by imposing that (10) be normaUzed to N . Thus is 
just the inverse of the level density Qnilj)'- 



1 



■ g4t,) = Jd^rg^iif). (12) 



The other quantity entering in the density of the Bose 
condensate, Eq. (10), is the chemical potential fj, which 
corresponds to the lowest eigenvalue of the GPE. In order 
to determine the chemical potential fi, a requantization 
of the TF approximation is necessary [15]. The need for a 
requantization of the TF theory for individual states has 
been recognized in Ref. [27] and our procedure of requan- 
tization clearly follows what is proposed there. The stan- 
dard semiclassical quantization procedure is given by the 
Wentzel-Kramers-Brillouin (WKB) method. However, in 
order to have a more explicit formula, we apply here the 
simplified method described in Ref. [15], which becomes 
exact in the three dimensional HO case. Thus for the 
non-interacting case we fix the chemical potential to be 
equal to the GPE eigenvalue, Eq. (5). 

To proceed further it is useful to write the Bessel func- 
tion in Eq. (10) as a power series [28]: 



(-1)* /ajx'^+si 



1=0 ^ ' 



(13) 



Using this result, performing the remaining integral, 
and remembering the identity 
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11,12=0 



2k + 2j 
j 



(14) 



we obtain the following expression for the local spectral 
density: 



mko (r) 



£ .t'r!':£^'r.. »-m. 



27r2fi2 j^^ j!(2K + j)\{2K + 2j + 1) 



(15) 

For the non-interacting harmonic oscillator the integral 
in Eq. (12) can be evaluated analytically, with the result 
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- =—y 



c« ^ iK2« + j)!(2K + 2j + 1)(2k + 1j + 2) 

(16) 

Fig. 1 displays the square root of the normalized TF 
density (10) for n = \ along the coordinate for z = 0, 
where HO units have been used. In the same figure the 
quantal wave function which describes the k = 1 vortex 
state [see Eq. (4)] is also plotted. As it can be seen from 
Eq. (15), for ^ the semiclassical TF density goes to 
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FIG. 1: Square root of the normalized TF ill 0) density 
and quantal wave function (GPE) for a non-interacting Bose 
condensate with k = 1 as a function of for 2 = 0. The third 
curve corresponds to Eq. (20). Wave function and radius are 
given in HO units (a^jo^ and uho, respectively). 



zero in the same way as the quantum mechanical result, 
i.e., Puir) oc r^'^. At the classical turning point [V{r) = 
/i] the TF density goes to zero as oc [A;o(r)]^+2K_ xhus 
the turning point will be changed by the interaction only 
via the change of the chemical potential /U. 



The reader unfamiliar with the Thomas-Fermi ap- 
proach may be worried about the locally relatively strong 
deviations of the semiclassical density from its quantal 
(GPE) counterpart. In this respect it should be remem- 
bered that the TF densities must be considered in the 
sense of distributions [26] [see e.g. the step function in 
Eq. (10)] and, therefore, they only make sense when 
used under integrals to calculate expectation values of 
"slowly varying" operators. In fact, the H corrections 
to the density (not considered here) still deviate much 
more strongly from the true quantal densities, since they 
contain a square-root singularity at the classical turning 
point. Nonetheless these h corrections, when used under 
integrals, improve the results for cases where the gradi- 
ents of the potential are not too strong [29] . It has been 
found in the past [15, 26] that when used in this way the 
TF approach (eventually with inclusion of H corrections) 
can yield very accurate results for expectation values. 
The relation of the TF approach with h"^ corrections to 
the WKB approach has been discussed in Ref. [27]. 



The TF kinetic energy density can also be derived from 
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the spectral density matrix (9) as 

f dk, dk 



= [fi-V{R)]g^M- 



(17) 



Using Eqs. (10) and (17) it can easily be checked that 
the expectation values of the kinetic and potential ener- 
gies fulfill the virial theorem as it is expected in the TF 
approach [30]. Thus due to our choice of the chemical po- 
tential jj, (5), our semiclassical TF approximation to the 
vortex state in the non-interacting case exactly repro- 
duces the quantal expectation values of the kinetic and 
potential energies in spite of the aspect of the semiclas- 
sical density profile as compared with the quantal one. 

It is easy to see that with the HO potential the argu- 
ment of the Bessel function entering in the density cannot 
become large: 



< ko{f}\r\ < 



JJ_ 



(18) 



For example, in the case k = 1 the argument becomes at 
most 5/2, and already the first four terms of the expan- 
sion (13) give an accuracy better than 0.5%. In the case 
K = 0, the result of Ref. [15] is recovered if one takes 
only the first term of the expansion in the TF density, 
Eq. (10). 

For completeness we note that in the literature also a 
different approach for projecting the semiclassical den- 
sity matrix onto good angular momentum and can 
be found [31]. Repeating the steps described there for 
the projection onto good only (i.e., essentially using 
asymptotic expansions for the Bessel functions) one finds 
the following expression for the Wigner transform of the 
spectral density matrix: 



gi:{R,p) = hd{H'^'-,,)S{Lf-hK) 



(19) 



with H"' = f/{2m) + V{K) and U^^ =Rxp. From this 
formula the density is easily obtained by integration over 
p: 

The constant is determined by the normalization 
condition, which for the non-interacting case results in 
Ck = 2h^uj'^ /{fi— nfiw) = 4hu;/3. The density profile cor- 
responding to Eq. (20) is also shown in Fig. 1. From this 
figure it is evident that Eq. (20) makes sense only as a 
distribution for the calculation of expectation values and 
not for the calculation of local quantities like the density 
itself. However, since the density given by Eq. (20) does 
not at all depend on the shape of the potential (except 



for the determination of the turning points), this form 
seems difficult to be used for a self-consistent calculation 
in the interacting case. Let us note that again the virial 
theorem is fulfilled and expectation values of operators 
can be obtained very accurately [31, 32]. 



III. THE INTERACTING CASE 

Let us now disciiss the TF approximation to the quan- 
tal solution of the GPE (2). Of course the semiclassical 
formalism described in the previous section can still be 
applied provided that the potential V{f) in the interact- 
ing case is given by 



(21) 



As it was mentioned before, the TF density correspond- 
ing to the vortex state (10) depends on two independent 
constants to be determined: the normalization and 
the chemical potential /i. In the interacting case, and 
following the same strategy as in Ref. [15], we determine 
Ck and iJ. by imposing that the TF density be normalized 
to the number of particles N in the Bose condensate and 
that the integrated level density 



d^r0[ii-V{f)] 



k'oir) 



E 

3=0 



27r2 

(-1)-'' [fco(r)rj2«+2j 



j!(2« + j)\{2K + 2j + 1)(2k + 2j + 3) 



(22) 



become equal to that of the non- interacting HO, which 
for jj. = (3/2 -I- K)hLj is given by 



N. 



HO 



= E 



{-ly {2k + 2j)\ (3/2 + 

j\{2K + jy.{2K + 2j + 3y. 



(23) 



The strategy for the self-consistent solution in the in- 
teracting case is now very simple. Instead of starting with 
a fixed particle number N, it is convenient to choose some 
value for TVc^. Then we choose some initial value for ^. 
For given Nc,^ and /x we solve Eq. (15) for g^{r), which 
is non-linear since also the right-hand side depends on 
gj^{f) through 



fco(r) = 



h 



(24) 



Then the integral (22) is evaluated and the result is 
compared with the corresponding result for the non- 
interacting harmonic oscillator, Eq. (23). If the level 
number is too small {N^ < iV^*^), /i is increased, oth- 
erwise {N^ > N^'-') IJ. is decreased. This procedure is 
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iterated until Nj^ = Nj^'-' . Finally the particle number is 
obtained by evaluating the integral 

N = Jd^r p{f) = Nc^ J d\ g^:{f) . (25) 

Before comparing the results obtained within our TF 
approach to the results from solving the GPE numeri- 
cally, let us briefly discuss two approximation methods 
which have been developed for the case of large iV. The 

first one, known as the TF limit in the literature and 
discussed, e.g., in Refs. [10, 11, 16, 17], is obtained by 
dropping the kinetic energy part Ckm coming from the 
radial and axial motion and retaining only the rotational 
part Brot of the total kinetic energy, i.e., only derivatives 
with respect to (p in Eq. (2). Under this assumption it is 
easily obtained that the N oo limit of the density of 
the vortex state reads 



(26) 



In this limit the density vanishes inside of r^min and out- 
side of rj^max; defined by the zeros of Eq. (26). The 
chemical potential fi is obtained through the particle- 
number condition. The formula (26) has the advantage 
that it represents an analytic expression for p{r), but it 
is clear that the result p = inside a vortex core with 
radius rj^min is not realistic. We will call hitherto formula 
(26) the iV ^ oo TF limit. 

The second approximation method, known as the 
"method of matched asymptotics" (MA), was introduced 
in Ref. [33] to describe the dynamics of vortices, and used 
in Ref. [34] to calculate the energy of a static vortex. We 
will follow here the simplified derivation for the case of 
a straight vortex given in Ref. [18]. First let us briefly 
review the description of a vortex state in a system with 
VexL = 0. In this case it is useful to define the asymptotic 
density po = p/g and the healing length = h/ y/2mpog 
[10], and to write the condensate wave function in the 
form 



</.(rJ = 




p ( ^/2mp 
g \ h 



(27) 



Inserting this expression into the GPE (2) with Vext = 
0, one obtains the following differential equation for the 
function /„: 

~fLi^) - a^) + ^u^) + fi{^) = u^) ■ (28) 

With the boundary conditions (for k > 1) 

/«(0) = and lim/,(x) = l (29) 

X — >oo 

this differential equation can be solved numerically [10, 
35]. Now we turn to the case of a trapped system. In 
the limit of large N it is is clear that the external poten- 
tial Vext can be regarded as constant on the length scale 



^0 corresponding to the size of the vortex core. (More 
precisely, the condition which has to be fulfilled reads 
Na/ano > 1, as it is the case for the N oo TF ap- 
proach.) Thus we obtain an approximate description of 
the trapped system by replacing the chemical potential p 
by a local chemical potential p — Vexti^)- Inside the clas- 
sically allowed region [Vea;t(^) < fj] the order parameter 
then takes the form 

m = ^l^^^f.{^^^^^^^r.). (30) 

As before, the chemical potential p is determined by the 
particle-number condition. 

Note that in the region far away from the vortex core, 
i.e., for 7^ » ^0, one can expand /k(x) in powers of 1/x. 
Using Eqs. (28) and (29), one obtains 



Ux) = 1 - 



2a;2 



8x^ 



(31) 



Inserting this into Eq. (30) one immediately recovers Eq. 
(26). However, this shows that Eq. (26) is not valid for 
^ ^ ^0; i-e., inside the vortex core. Another difference 
between Eq. (30) and Eq. (26) concerns the behavior 
of the wave function at the outer classical turning point. 
In contrast to the usual N ^ oo TF limit, the kinetic 
energy corresponding to the wave function (30) is not di- 
verging, since the square-root — Vext in Eq. (30) is 
multiplied by the function /„, which is proportional to 
{p — VextY^"^ near the classical turning point. Neverthe- 
less it is not reasonable to use Eq. (30) to calculate the 
kinetic energy near the turning point, since the decrease 
of the function i s just indicati ng that the local healing 
length ^(r) = h/y/2m{p — Vext) becomes large and that 
the approximation breaks down. 

Now we proceed to a detailed numerical comparison of 
the (ft 0) TF predictions with the exact quantal values 
obtained from the GPE, Eq. (2). For our numerical ap- 
plication we consider ^''Rb atoms in a spherical trap rep- 
resented by a HO potential with length ano = 0.791 pm 
[9]. The s-wave scattering length is taken as a = 100 ao 
[6] where ao is the Bohr radius. Table I collects the 
chemical potential (/i), the total {ctot), HO (ej^o), self- 
interaction (csei/), and kinetic energies per particle for 
vortex states of condensates with 100, 10^, and 10^ atoms 
in the trap. The kinetic energy is split into the rotational 
part {crot) and in the one corresponding to the radial and 
axial motion {ekm)- The numerical values displayed in 
Table I show that our (ft 0) TF approach reproduces 
very well the quantal eigenvalue {p) as well as the total 
energy per particle {etot) even for a small number of par- 
ticles such as 100. The agreement between the quantal 
and TF values improves when the number of particles in 
the condensate increases, as it is expected. The HO and 
the self-interaction contributions to the total energy are 
also well reproduced by our semiclassical approach. For 
very large numbers of particles (TV = 10^) the quantal 
results are also well reproduced by the N ^ oo TF limit. 
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N 



etot &HO eself Crot 



GPE 


2.74 


2.62 


1.34 


0.12 


0.480 


0.686 


N ^ oo 


2.72 


2.61 


1.33 


0.10 


0.501 


0.673 


1.88 


1.59 


0.87 


0.29 


0.438 





MA 


1.86 




0.84 


0.23 


0.689 





GPE 


8.40 


6.30 


3.67 


2.10 0.271 


0.255 


10^ - 
N ^ oo 


8.28 


6.19 


3.62 


2.09 0.350 


0.130 


8.19 


5.99 


3.54 


2.19 0.253 





MA 


8.23 




3.57 


2.17 0.272 





GPE 50.18 35.93 21.53 14.26 0.087 0.059 

h ^O 50.13 35.86 21.50 14.27 0.116 -0.024 

N ^oo 50.14 35.86 21.50 14.28 0.083 

MA 50.14 — 21.50 14.27 0.086 — 



10' 



TABLE I: Chemical potential (/u) and energy per particle 

(etot) and its different contributions in hcu units: harmonic 
oscillator energy (cho), interaction energy (6,^// ), and kinetic 
energy split into its rotational (crot) and radial and axial 
(efcin) parts. The parameters chosen correspond to a single- 
quantized vortex (k = 1) in a spherical trap [ano = 0.791 /nm) 
containing 100, 10*, and 10^ ^'^Rb atoms (scattering length 
a = 100 oo). The results obtained from the GPE are com- 
pared with the results from the (ft 0) TF approach and 
from two approximation methods for large A^: the so-called 
AT — »■ oo TF method, Eq. (26), and the method of matched 
asymptotics (MA), Eq. (30). Note that eun is neglected in 
the N oo TF limit and not accessible within the matched- 
asymptotic approach. 



Eq. (26), because the neglected contribution (i.e., the 
kinetic energy due to the radial and axial motion Ckin) 
is very small. However, it should be pointed out that 
the key assumption of this N ^ oo limit is not fulfilled, 
because the kinetic energy of the radial and axial mo- 
tion is still of the same order as the rotational energy, 
as can be seen from the quantal results (GPE) listed 
in Table I. In fact, even in the limit N ^ oo the ra- 
tio ekin/^rot does not go to zero (see appendix). The 
method of matched asymptotics, Eq. (30), gives better 
results than the N ^ oo TF limit, Eq. (26), except in 
the case of small numbers of particles {N = 100), where 
both large- methods fail. 

Concerning the kinetic energy some comments are in 
order. First of all, we want to point out that the N ^ oo 
theory neglects (and in fact cannot access [15]) the con- 
tributions coming from the radial and axial motion, so 
they are not listed in Table I. For a small number of 
atoms, such as 100, our limit is able to reproduce 
reasonably well both, Crot and ekin contributions to the 
total kinetic energy per particle. When the number of 
the atoms in the trap grows, the total kinetic energy per 
particle decreases and the agreement between the quan- 
tal result and the TF prediction worsens for this quan- 
tity. This situation is also found in the ground-state case 
discussed in Ref. [15] where the (small) quantal and TF 
kinetic energies can differ by a factor two for a large num- 
ber of particles (see Table II of Ref. [15]). The reason for 



these disagreements between the quantal and TF kinetic 
energies for large number of atoms in the condensate lies 
in the fact that in this case the kinetic energy is domi- 
nated by quantal corrections that are non-analytical in h 
and consequently cannot be reproduced in a pure TF ap- 
proximation [15]. A detailed comparison shows that the 
{h 0) TF theory systematically overestimates the rota- 
tional part and underestimates the axial and radial parts 
of the kinetic energy, the latter even becoming negative 
for very large numbers of particles, although the total 
kinetic energy remains always positive. The reason for 
this behavior is that the TF density is too high inside 
the vortex core, as will be discussed below. 

It should be pointed out that, as happens for the non- 
interacting case, the virial theorem, which for the inter- 
acting case reads [6] 



2(efci„ -I- Crot) — '^e.HO + 3esei/ — , 



(32) 



is also fulfilled in our TF approach to vortex states for a 
Bose condensate in a spherical trap. 

Figures 2-4 display the normalized order parameter for 
100, 10"*, and 10^ atoms of *'^Rb in the trap along the 
radial axis for 2; = 0. The dashed line corresponds 
to the {h 0) TF limit. For comparison we show the 
corresponding order parameter obtained from the quan- 
tal solution of the GPE (2) (solid line), which is obtained 
through imaginary time step techniques [6] , and the order 
parameter obtained from the method of matched asymp- 
totics, Eq. (30) (dashed-dotted lines). Looking at the 
shape of the semiclassical {h 0) compared with the 
quantal order parameter one can see that the agreement 
increases with the number of particles in the condensate, 
as it happens in the TF approximation for the ground 
state. The effect of the self-interaction that progressively 
modifies the density profile of the condensate in the vor- 
tex state with respect to the non-interacting case is also 
followed by our semiclassical TF densities. Only inside 
the vortex core (r^ « 0) the agreement worsens with in- 
creasing number of particles. 

This can easily be understood by looking at the corre- 
sponding self- consistent potentials shown in Fig. 5. The 
main assumption of our semiclassical TF theory is that 
gradients of the potential can be neglected. This assump- 
tion becomes more and more justified with increasing 
number of particles, except in the vicinity of the z axis 
{r^ ~ 0), where the self consistent potential rises rapidly 
from zero to w fi. For the case of moderate numbers of 
particles, the semiclassical description of the vortex core 
could be improved by considering higher h corrections to 
the TF solution, which take into account the gradients of 
the potential. However, we should remember that the h 
or gradient expansion is an asymptotic series which can 
only work as long as the gradients of the potential are 
not too strong, even though the theory often works quite 
far beyond its limits (see Ref. [29]). For very steep po- 
tentials only a partial resummation of the h series like 
in WKB, to account for the nonanalytical behavior in h, 
can help. This will further be discussed in the Appendix. 
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FIG. 2: Normalized order parameter obtained from the GPE, 
from the (h —> 0) TF limit, and from the approximation of 
matched asymptotics (MA) for large [Eq. (30)], of an in- 
teracting Bose condensate of 100 ^'^Rb atoms in a spherical 
trap with ano = 0.791 fim in a vortex state with k = 1 as a 
function of for 2: = in HO units. 



FIG. 5: Self-consistent potentials in units of hcu as a function 
of rj_ in units of uho obtained in our {h — > 0) TF approach 
corresponding to the density profiles shown in Figs. 2-4. 
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FIG. 3: Same as Fig. 2, but for 10* atoms in the trap. 
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FIG. 4: Same as Fig. 2, but for 10*^ atoms in the trap. 



As can be seen in Figs. 3 and 4, for large numbers 
of particles the density profiles obtained from Eq. (30) 
follow remarkably well the quantal profile except near the 
classical turning point where the approach breaks down. 

Comparing the so-called N 00 TF approach [Eq. 
26)] with the h ^ TF approach proposed in this work, 
one can see from Table I that our TF method reproduces 
better the different quantal contributions to the energy 
of the vortex state for small (A^ = 100) and moderate 
{N = 10'*) numbers of particles in the condensate while, 
for large numbers both limits (^ — > and N — > 00) co- 
incide. Concerning the density profiles, the difference is 
obvious: In our approach the density profile goes like 
(X as in the quantal case, whereas in the N ^ 00 
TF limit there is a small p = region determined by the 
inner turning point r^min- It should be mentioned that 
within the improved (matched asymptotics) N ^ 00 ap- 
proximation [Eq. (30)], the density profile also goes like 
y/p cx r^, and contrary to our ?l — > approach it is ca- 
pable to reproduce the density profile in the vortex core 
in the case of large A^. However, this method has noth- 
ing to do with the semiclassical asymptotic h expansion 
considered here. 

Finally it should also be pointed out that our TF limit 
is able to deal with vortices in the attractive case (neg- 
ative scattering length). In this case the kinetic energy 
is crucial and the large- A^ limit is not well-defined. The 
same is true for the description of the ground state (i.e., 
no vortex), as shown in Ref. [15]. There the h 
approach has been used in the repulsive as well as in 
the attractive case, whereas the N 00 approximations 
(the so-called AT — > 00 TF limit as well as the method of 
matched asymptotics) can be applied only in the repul- 
sive case. 
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IV. APPLICATION TO VORTICES IN 
SUPERFLUID TRAPPED FERMIONIC GASES 

In this section wc will describe how our TF approach 
can also be used for the description of vortices in super- 
fluid fermionic systems. This is possible since at least for 
a certain range of temperatures, the so-called Ginzburg- 
Landau regime, the order parameter A(r) is described by 
an equation which has exactly the same form as the GPE 
(2). As derived in Ref. [24], for temperatures T near the 
critical temperature and for low trapping frequencies 
uj {hu> <C ksTc) the Ginzburg-Landau equation (GLE) 
reads 



1 + 2A 



.(0) 



2A 



+ 



p2 

7C(3) 
87r2 



In- 



A(r) 



A(f) 



A(r)=0, (33) 



with the definitions K = ^7C{3) / {48tt^) huj/{kBT), 
A = 2fcF|a|/7r, and Rtf = hkF/imu!), where kp de- 
notes the local Fermi momentum at the center of the 
trap. The temperature Ti°^ is the critical temperature 

of a homogeneous system having the same density as 
the trapped system has at the center. It is given by 



n{0) 



(8e 27/7r)e_Fe 



-l/A 



[36], with 7 « 1.781 and 



ep = h'^kj,/{2m). 

It is convenient to rewrite Eq. (33) in terms of dimen- 
sionless quantities. To that end we define 



1 \l/2 



1 



1 \i/4 r* 



KJ V 2XJ Rtf 
7C(3) 1 / 2A \i/2 



167r2 ii: Vl-f 2A/ 
. 1 / 2A \i/2, Ti°^ 



keT 



(34) 
(35) 

(36) 
(37) 



With these definitions, Eq. (33) becomes 



+ )f + = fim , (38) 

which is the same as the GPE rewritten in HO units, i.e., 
with the replacements f/uHO r, g/{hLoa%Q) g, and 

4"^ HO ^ V- 

However, there is one important difference between the 
GPE describing the Bose-Einstein condensate and the 
GLE describing the order parameter A(r) of a superfluid 
Fermi system. In a Bose-Einstein condensate, the par- 
ticle number N , i.e., the norm of is fixed, and the 
chemical potential jl has to be determined from the GPE 
(38). For the GLE the situation is reversed: The chem- 
ical potential fi is fixed by the temperature T and other 
parameters [Eq. (36)], whereas the normalization of 
i.e., the magnitude of the gap A, has to be determined 
fromEq. (38). 



The lowest possible value of jl, for which a solution 
of the GLE (38) can be found, corresponds to the case 
that the normalization A'' goes to zero, such that the 
non-linear term can be neglected. In this case Eq. 
(38) reduces to the Schrodinger equation of the non- 
interacting harmonic oscillator with the lowest eigenvalue 
Amin = 3/2. This gives an upper limit for T/TjP\ which 
was used in Ref. [24] to estimate the critical temperature 
Tc of the trapped Fermi system. 

In this article we arc interested in vortex states, i.e., 
in solutions of the form (1). In the framework of the 
GL theory, vortex states of superfluid Fermi systems are 
discussed in Ref. [25], where the GLE (38) is solved for 
a two-dimensional geometry, i.e., (/)(r^,z) = (/>(?^), cor- 
responding to a trap with an extremely elongated po- 
tential. In two dimensions the lowest possible value of 
jl, for which vortex solutions can be found, is given by 
Mmirt ~ 1 + However, as in our discussion of vor- 
tex states in Bose-Einstein condensates, we will con- 
sider the spherical case, in which the maximum temper- 
ature for the existence of vortex states is determined by 

jj"mm — 3/2 -|- K. 

Since the GLE is identical with the GPE, it is obvi- 
ous that our TF approach described in the previous sec- 
tions can immediately be applied also to the GLE. Only 
the iteration procedure for the self-consistent solution is 
somewhat different, since now ji is given instead of N . 
We start with some guess for TVc^ and calculate the in- 
tegrated level density Eq. (22). Now, if N'^ < N^'^ , 
the value of Nck is increased, otherwise it is decreased. 
This procedure is iterated until A^^ ~ ■ Due to this 
quantization rule it is clear that A goes to zero when the 
temperature approaches the critical temperature corre- 
sponding to jl = 3/2 -|- K. 

We are now going to compare the results of our TF 
approach with the fully quantal solution of Eq. (38). 
The parameters used for our calculations are taken from 
Ref. [22], i.e., wc consider A^en = 573000 ^Li atoms 
(scattering length a = — 2160ao) in a trap with a; = 
27r X 144 Hz. The self-consistent mean-field potential 
of the cloud has been neglected in the derivation of the 
GLE (33) in Ref. [24], but we take it into account in 
an approximate way by replacing the external trapping 
frequency a; by a higher one, ujejj = 27r x 170 Hz, as 
it has been done, e.g., in Ref. [23]. The parameters 

Rtf and Ti°^ are obtained from ep = {3N6i^i)^/^hu!eff 
and the relations given below Eq. (33), with the result 

Rtf = 48.7 /xm and Ti°' = 36.7 nK. The temperature 
corresponding to /Et = 3/2, i.e., the critical temperature 

of the trapped system, is Tc = 31.2 nK. 

In Fig. 6 we show the order parameter A obtained 
from the numerical solution of Eq. (33) (sohd lines). For 
the parameters listed above, the lowest temperature for 
which vortex states can exist, i.e., the temperature for 
which jl = 5/2, is approximately 0.86 T^. Therefore we 
display the order parameter only for temperatures below 
this value, namely for T/Te = 0.85, 0.8, and 0.75. For 
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FIG. 6: Order parameter A of a supcrfluid trapped Fermi gas 
in a vortex state with /t = 1 as a function of for 2 = 0. The 
parameters were chosen corresponding to 573000 ^Li atoms in 
a spherical trap with a; = 27r x 144 Hz. The order parameter 
A is given in units of fesTc (Tc = 31.2 nK), the radius in 
units of Rtf {Rtf ~ 48.7 ^m). SoUd hnes result from the 
numerical solution of the GLE (33), whereas the dashed linos 
are obtained within the (/t 0) TF approach. The three pairs 
of curves correspond to three different temperatures (from 
bottom to top: 0.85Tc, 0.8rc, and 0.75Tc) 



T/Tc 




N 


Fgl {pK) 


0.85 


GLE 


2.12 
2.37 


-0.0079 
-0.0088 


0.8 


GLE 


16.37 

18.41 


-0.391 

-0.440 


0.75 


GLE 34.18 
38.20 


-1.44 
-1.61 



TABLE II: Normalization N of the order parameter and GL 
free energy Fql for the parameters used in Fig. 6, obtained 
from the numerical solution of the GLE and from our {h — > 0) 
TF approximation. 



T/Tc = 0.85 the order parameter is still very small, and 
it increases rapidly as the temperature decreases. As can 
be seen in Fig. 6, in all three cases the "amplitude" 
of the order parameter and the position of the maximum 
are well reproduced by our TF approximation. Note that 
also the vortex core is well described. 

As already stated, our TF solution has to be inter- 
preted in terms of distributions, and then the agreement 
is even better than it seems from Fig. 6 if one looks 
at integrated quantities. As an example we consider the 
normalization N = J (fif\(j){f^\^. For the three tempera- 
tures mentioned above, the normalizations obtained from 
the numerical solution of the GLE and those correspond- 
ing to our TF approximation are in good agreement, as 
shown in Table II. As a more meaningful example for 
an integrated quantity let us look at the GL free energy 
Fgl- The explicit expression for the functional Fgl[^] 



is given in Ref. [24]. Following this reference, we retain 
only the leading terms in the small quantities K, r/ Rtf, 
ln(Ti°Vr), and A/{kBTc). Then, after integration by 
parts, the GL free energy functional becomes 

Fgl[A] = ^ld'r{- K-Rl^A*V-A 

7C(3) 



+ 



703) ^_ 4-1 .39. 



167r2 {UbT) 

In the TF approach, the first term (oc A* V^A) cannot be 
obtained directly from the TF approximation for A(r), 
but it rather has to be calculated analogous to the kinetic 
energy density in Eq. (17). As a consequence, most of the 
terms cancel, as it is the case if A(r) is the exact solution 
of the GLE (33), and only the last term (oc | A|^) survives, 
but with negative sign. Thus, for the TF approximation 
as well as for the exact solution of the GLE, we can write 
in terms of the dimensionless variables defined above 



Fgl = 



4e|(fcBT)2 



K5/2 



(l 



2A \i/4 



-I-2A 



(-fw). 



(40) 

Results for Fgl obtained from the numerical solution of 
Eq. (33) and from the TF solution are listed in Table 
II. The agreement is as good as for the normalizations 
TV. In fact, the deviations are mainly due to the different 
normalizations, i.e., the ratio Fgl/N obtained from the 
TF approximation is very close to the exact one. 

To conclude this section, we stress that in the range 
of validity of the GLE the "chemical potential" jl must 
not become large. This is the reason for the rather small 
normalizations of the order parameter and results in a 
shape of the order parameter as a function of r which 
resembles very much the shape of a non-interacting HO 
wave function. Under these conditions it is clear that the 
N ^ 00 limit cannot be used as an approximate solution 
of the GLE, as has also been noted in Ref. [25]. 



V. SUMMARY AND CONCLUSIONS 

Using the semiclassical Thomas-Fermi approximation 
understood as h ^ limit rather than N ^ 00 limit, we 
have studied the vortex states of a Bose condensate of 
atoms confined in a spherical magnetic trap. 

We started analyzing the vortex states in a non- 
interacting trapped Bose gas. Due to the symmetry of 
the problem, we have obtained first the Thomas-Fermi 
density projected on states of defined i^. In this non- 
interacting case the density is normalized by adjusting 
the normalization constant c^, and the chemical poten- 
tial /i is fixed, according to the WKB quantization rule, 
to the quantal eigenvalue of the quantum state. 
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In the interacting case the norniahzation constant and 
the chemical potential arc fixed to normalize the Thomas- 
Fermi density to the number of particles and that the 
integrated level density become equal to that of the non- 
interacting case. For particle numbers where the kinetic 
energy coming from the radial and axial motion is a non- 
negligible part of the total kinetic energy, our Thomas- 
Fermi approach, understood as h ^ limit, yields very 
satisfying results as compared with the corresponding 
quantal values. For a very large number of particles in 
the condensate, the small Thomas-Fermi kinetic energy, 
obtained in our approach, is smaller than the quantal ki- 
netic energy, which, for large number of particles, is also 
dominated by quantal corrections as it happens for the 
ground state [15]. 

The vortex state density profiles obtained in our 
Thomas-Fermi approximation reproduce quite well the 
quantal ones, especially for a very large number of parti- 
cles. However, inside the vortex core our Thomas-Fermi 
densities are too high. Also near the classical turning 
point our TF densities locally fail because at this point 
the density is completely dominated by quantal contri- 
butions which are non-analytical in h and which cannot 
be reproduced by semiclassical approximations of the TF 
type. However, it shall be kept in mind that the semi- 
classical density has to be understood as distribution very 
efficient for describing expectation values rather than lo- 
cal quantities such as the density profile. In this sense 
we see that the quantities presented in Table I are much 
more accurate than one would expect from an inspection 
of the local densities shown in Figs. 2-4. 

The approach is also well suited for the description 
of vortex states of superfluid trapped fermionic systems 
in the GL regime, where the various approximations de- 
veloped for large N cannot be used at all. It should be 
mentioned that the conditions for the validity of the GLE 
imply that the parameters of the equivalent GPE always 
correspond to a rather small number of particles. In this 
case the normalization of the order parameter (see Table 
II) and the position of the maximum are well reproduced 
by the — > limit as compared with numerical solutions 
of the GLE. Also the vortex-core region is well described 
in this case. 



APPENDIX: LARGE- iV LIMIT FOR VORTEX 

STATES 

In order to discuss the large- iV limit for vortex states 
more thoroughly, we start from Eq. (30), which, as shown 
in Fig. 4, becomes very accurate in the limit of large N. 
Let us look at the different contributions to the kinetic 
energy, Crot and etin- In the infinite system, the ener- 
gies per unit length, dErot/dz and dEkm/dz can easily 
be obtained from the numerical solution for /k. Since 
dErot/dz diverges logarithmically, the corresponding in- 
tegral has to be cut off at some radius R [35]. For k = 1 
the results read (i? » ^o): 



dEr. 



dz 
dEki, 



' ' In — - 0.40 



dz 



= 0.28 



m \ ^0 



m 



(A.1) 
(A.2) 



In complete analogy to the derivation of the total energy 
of a vortex in a trapped system in Ref. [37], one can 
use these results to obtain explicit expressions for the 

rotational and radial kinetic energies of a vortex, Crot 
and 6^°^*^, which for a spherical trapping potential read 



N 6 m, \ 



In- 



6 



1.18 



-'kin 



28-^-r 
N 3 m 



(A.3) 
(A.4) 



Here Vmax = -v/2/i/(mw^) is the radius of the condensate. 

However, the kinetic energy has also another contribution 
due to the finite size of the trapped system. Since 



trap 
''kin 



outside the vortex core the shape of the condensate is 
almost not changed, we assume that for this contribution 
the relation derived in Ref. [38] for the case without 
vortex, remains valid: 



-'kin 



2 mrj, , 



In- 



).26) . 



(A.5) 



Since the volume of the vortex core is negligible in the 
limit N ^ CO, j2 depends on N in the same way as in the 
large- A'' limit for the ground state, 
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1^ 



hiJ/15Na\^/5 



2 V aHO ' 
Using this, we finally obtain 

e™t = ?iw In 2.95 



^kin 



a-HO 

\lhNa) \2 aHO ' 



(A.6) 

(A.7) 
(A.8) 



From these equations we conclude that the ratio 
^kin/^rot does not go to zero, but approaches 1/2 for 
N oo. Hence, neglecting the radial and axial parts of 
the kinetic energy, but retaining the rotational part, as 
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it is done in the literature [10, 11, 16, 17], is not justified 
and does not correspond to the proper iV ^ oo limit. 
Instead, the correct large- limit is given by Eq. (30), 
except at the surface of the condensate. The latter can 
be approximated, e.g., by the exact solution of the GPE 
for a linear potential, as it has been done in Ref. [38] in 
order to derive Eq. (A. 5), and also in Ref. [37]. 

It should, however, be noted that Eq. (30) and Eq. 



(A. 5) correspond to a partial resummation to all orders 
in h as demonstrates the nonanalytical dependence on h 
of these quantities. Such resummation techniques, also 
encountered in the WKB approximation, are necessary 
whenever the asymptotic Wigner-Kirkwood h expansion 
breaks down. This is always the case when the gradients 
of the potential start to diverge like in the vortex core 
for A'' —> 00, see Fig. 5. 
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